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Abstract 


Theories of modified gravity suggest that the propagation speed of gravitational waves (GW) v, may deviate from 
the speed of light c. A constraint can be placed on the difference between c and v, with a simple method that uses 
the arrival time delay between GW and electromagnetic wave simultaneously emitted from a burst event. We 
simulated the joint observation of GW and short gamma-ray burst signals from binary neutron star merger events in 
different observation campaigns, involving advanced LIGO (aLIGO) in design sensitivity and Einstein Telescope 
(ET) joint-detected with Fermi/GBM. As a result, the relative precision of constraint on v, can reach ^10 7 
(aLIGO) and ~10 !5 (ET), which are one and two orders of magnitude better than that from GW170817, 
respectively. We continue to obtain the bound of graviton mass m, < 7.1(3.2) x 10 ??eV with aLIGO (ET). 
Applying the Standard-Model Extension test framework, the constraint on v, allows us to study the Lorentz 
violation in the nondispersive, nonbirefringent limit of the gravitational sector. We obtain the constraints of the 


dimensionless isotropic coefficients s(? at mass dimension d = 4, which are —1 x 107^ < s? <9 x 10-7 


for aLIGO and —4 x 10716 < 3$ < 8 x 


10-3 for ET. 
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1. Introduction 


Gravitational waves (GWs), which are perturbations created 
by the bending of spacetime due to the distribution of mass, are 
one of the key predictions of Einstein's theory of General 
Relativity (GR). The first direct detection of a GW event, 
GW150914 (Abbott et al. 2016), was observed from a binary 
black hole (BBH) merger during the first observing run (O1) of 
Advanced Laser Interferometer Gravitational-Wave Observa- 
tory (aLIGO) (Aasi et al. 2015), opening a new window to 
astronomy and physics. The GW observations enable us to test 
gravity theories in strong and dynamic regimes of gravity. To 
test gravity theories with GWs, it is crucial to search for 
anomalous deviation from the prediction of GR in model- 
independent ways. One of the tests is measuring the 
propagation speed of a GW. The physical implications of 
limiting v, are manifold. The limit on the speed of GWs is 
closely related to the validity and accuracy of general relativity. 

According to GR, in the limit in which the wavelength of 
GW is small compared to the radius of curvature of the 
background spacetime, the waves propagate along null 
geodesics of the background spacetime (Will 2014), i.e., a 
GW propagates with the speed of light c in vacuum. In the 
theories of modified gravity, v, could deviate from c due to the 
modification of gravity such as the coupling of gravitation to 


"background" gravitational fields. A GW generically propa- 
gates at a speed different from c and with frequency 
dependence dispersion relations in some theories of modified 
gravity (see, Kostelecky & Mewes 2016; Blas et al. 2016; 
Yunes et al. 2016; Bettoni et al. 2017; de Rham et al. 2017). On 
the observation with the GW detector, the propagation final 
resulted in an offset in the relative arrival times of 
simultaneously emitted messengers (GW and neutrinos, or 
photons, etc.) at the Earth (Will 2014; Lombriser & 
Taylor 2016; Abbott et al. 2017b). 

Currently, the most exact way to measure v, is by measuring 
the time delay between GW and electromagnetic (EM) signals 
of the same astrophysical source (Cornish et al. 2017). With the 
Advanced Virgo (Acernese et al. 2014) detector joining in later 
O2 (Abbott et al. 20192), the coalescence of a binary neutron 
star (BNS), GWI170817 (Abbott et al. 2017b, 2017c) was 
observed by the three-detector network (Abbott et al. 2020), 
which is a significant observational discovery in the field of 
multi-messenger astronomy. In Abbott et al. (2017b), LIGO 
and Virgo Scientific Collaboration (LVC) has constrained the 
difference between v, and c to be between (—3 x 107 9c, 
4-7 x 10 6c) by using the distance to the astrophysical source 
(the identification of NGC 4993 as the host galaxy of 
GWI170817 and GRBI170817A (Abbott et al. 2017a; 
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Coulter et al. 2017a, 2017b)) and the observed time delay of 
(+1.74 + 0.05) s between GRB170817A and GW170817. The 
observation of GW170817 also takes significant scientific 
discoveries that have set stringent limits on fundamental 
physics (Baker et al. 2017; Abbott et al. 2019b) and the origin 
of short Gamma-ray Bursts (more described in detail about the 
observations of GW170817 and GRB 1708174, see (Abbott 
et al. 2017c; Goldstein et al. 2017)), demonstrating the 
importance of multi-messenger astronomy. 

The LIGO-Virgo-KAGRA (LVK) 's O4 has been officially 
launched and has been running for the last few months. More 
scientific merits in fundamental physics, astrophysics, and 
cosmology can be excavated from a population of joint GW/ 
EM counterparts detection. It also facilitates the more stringent 
constraints on v,. According to scientific observation runs LVC 
O1-03 (Abbott et al. 2019a, 2023, 2024) and simulative studies 
(such as (Hendriks et al. 2023)), such a joint detection of GW/ 
GRB are rare (one out of ~10 binary neutron star (BNS) 
merger events). In this paper, we present a simulation study on 
constraining v,, using GWToolbox? (Yi et al. 2022; Hendriks 
et al. 2023) to simulate the joint GW-GRB detection of a BNS 
merger event with both second and third-generation ground- 
based GW detectors, including aLIGO in design sensitivity 
(Aasi et al. 2015) (design aLIGO) and Einstein Telescope (ET) 
(Sathyaprakash et al. 2010). In comparison to similar 
simulative work in Nishizawa (2016), we simulated the 
parameterized BNS population, which is consistent with the 
GWTC-3 (Abbott et al. 2023) catalog. 

We continue to place constraints on Lorentz invariance 
violation. A theoretical framework known as the Standard 
Model extension (SME) provides a comprehensive test frame- 
work for constraints on local Lorentz violation. The coefficients 
for Lorentz violation in the gravitational sector of the SME 
cause the modification of the group velocity of GWs. Some 
early work can be seen (Kostelecky & Mewes 2016; Colladay 
& Kostelecky 1998; Kostelecky 2004; Bailey & Kost- 
elecky 2006; Haegel et al. 2023). The constraint of v, allows 
constraining the nine coefficients for Lorentz violation in the 
limit of nonbirefringent, nondispersive at mass dimension 
d — 4. Thus, we use this bound to place a constraint on local 
Lorentz invariance violation within the context of the effective- 
field-theory-based test framework provided by the SME. 
Another way in which the speed of GWs could differ from c 
is if gravitation were propagated by a massive field (a massive 
graviton). The massive particle has a dispersion relation, 
FE? = m; + p? (in natural unit). Therefore, the bound on v, can 
be converted to a bound on the mass of graviton m, (Nishizawa 
& Nakamura 2014). 


5 Gravitational Wave Universe Toolbox: a Python package to simulate 
gravitational wave detections with different backgrounds. The web application: 
https:/ /gw-universe.org/, Python module: https:/ /bitbucket.org /radboudradiolab / 
gwtoolbox/src/master/ 
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The structure of this paper is as follows. In Section 2, we 
discuss the methods used to calculate the bound of v,, and 
Section 3 presents the constrain results of GW speed. In 
Section 4 and Section 5 we utilize the result of the GW speed to 
constrain both the coefficients for Lorentz violation and the 
graviton mass. 


2. Method 


In this section, we will discuss the methodologies employed 
to obtain constraints on vg, as well as our approach to 
simulating the observation on a BNS merger population using 
the GWToolbox. 


2.1. Propagation Speed of GW 


In the following, we are going to present a general method to 
constrain vg. We will start with a brief introduction to the 
method, comparing the arrival times between GW and a high- 
energy photon from a short Gamma-ray burst (SGRB) (Abbott 
et al. 2017b) to measure v,. 

Considering the compact binaries (like BNS) within the 
cosmological distances, assuming there is a small discrepancy 
ót between the propagation time of GWs and EM waves 
emitted by the sources, after traversing a distance of D. The 
fractional speed difference between GWs and EM waves (or 
photons) during the trip can be simply expressed as (Abbott 
et al. 2017b) 


— x c—, (1) 


We define óc = vg—c, is the difference between v, and c, the 
time decay ótz Afgy—-Af;,, where Atop, and Atin are the 
differences we observed in the arrival time and the intrinsic 
emission time difference of the two signals, respectively. 
Assume that the refractive index of the vacuum is unity, and 
both waves are expected to be affected by background 
gravitational potentials in the same way. This relation exhibits 
better constraint when considering sources at larger distances. 
However, for sources at cosmological distances (z > 0.1), the 
impact of cosmological effects cannot be disregarded. Third- 
generation ground-based GW detector such as ET enables us to 
observe NS-NS binaries at cosmological distances up to z ^ 2, 
and up to z~4 for NS-BH binaries (Sathyaprakash et al. 
2010). 

Now we extend the scenario to sources at cosmological 
distances, taking into account the effects of cosmological 
redshift and the influence of a specific cosmological model. 
Here we refer to the approach in Nishizawa & Nakamura 
(2014) and Nishizawa (2016) that considers the case of 
propagation over cosmological distances. Consider a back- 
ground in a flat lambda cold dark matter (ACDM) universe for 
simplicity. The comoving distance from the Earth to a source at 
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redshift z can be written as 


I z Vg EN óc 
x(0, z) J, qo" = ( : Ja. a Q 


The Hubble parameter H(z) is given by 


A(z) = Hoy Q1 F zy +O, (3) 


where €), and Qa are the matter density parameter and the 
cosmological constant, respectively, and Ho is the Hubble constant. 
Here use the cosmological parameters of the cosmos from Planck18 
(Aghanim et al. 2020), which are Hg—67.66km Mpc !s |, 
€),,, = 0.30966. We express Xo(0, z) as the comoving distance when 
v, — c. Consider a scenario in which GW and EM signals are 
emitted from the same source of redshift z with velocities of c and 
Vp, respectively. Both signals travel the same comoving distance 
from the source, yet they reach us at redshifts z= 0 and z= — Az, 
respectively. Neglecting the second-order (and more) correction in 
Az (more detailed derivations are given in Nishizawa & Nakamura 
(2014) and Nishizawa (2016)), for x(0, z + Az) = xo(0, z) gives the 
time delay ôt between GWs and photons introduced by óc/c is 
(Nishizawa & Nakamura 2014) 


Az óc [* dz 
ót — E ; 4 
Ho c J, H(z) a 


In addition, after redshifting the intrinsic time delay, the time 
delay in propagation ót becomes: 


6t = Atos — (1 + z) Afin, (5) 


where z is the redshift of the source. 

Here we use Equation (4) to obtain the fractional speed 
difference óc/c for any sources that know the distance or 
redshift and the time delay between the GW signal and EM 
signal in propagation. While v, is time-dependent in general 
(Bellini & Sawicki 2014; Saltas et al. 2014), and may deviate in 
varying scales from the speed of photons through coupling with 
the different gravitational fields during propagation in the 
current epoch of the Universe (if there is any modification of 
gravity which allows for a change in the GW propagation 
speed), in our work, we regard the speed of the GW as constant 
during propagation for simplicity, or rather, v, being measured 
here represents the average speed during a trip. 


2.2. Simulating the Observation on a BNS Merger 
Population 


As in the observations (Abbott et al. 2017b) of GW170817 
and GRB70817A, in our simulations, it is assumed that the GW 
signal arrives at the Earth before the GRB signal. To obtain 
more data on joint events that were heard (GW signal) and seen 
(EM signal) later, and to strengthen further constraints on the 
bounds of óc/c, we use GWToolbox to simulate the joint 
detection of GW and GRB by designing aLIGO and ET on a 
BNS merger population, using the same cosmological model 
used in Equation (4). The GRB detector for simulating joint 
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Table 1 
One Realization Based on BNS Merger Events Joint-detected by Design 
aLIGO and Fermi/GBM Over a 10 yr Observation Period with S/N 
Threshold — 8.0 


z dz ml m2 Atobs 
1 0.11 0.05 1.76 1,79 1.04 
2 0.06 0.02 1.96 1.22 0.99 
3 0.08 0.04 1.60 1.61 1.00 
4 0.06 0.03 1.31 2.24 0.85 
5 0.07 0.03 1.42 2.47 1.25 


detection here utilizes only the Fermi/GBM. A^fto»s is explained 
to be the time difference between the peak (usually in the 
merger phase) of the GW signal and the first sGRB photon in 
the propagation. Using the joint observation module of the 
GWToolbox, after setting a campaign of detection, one can get 
the catalog that includes a list of parameters of the samples 
where the BNS merger is successfully joint-detected by both 
the specified GW and GRB detectors in all simulated GW event 
samples. This includes the masses of the binary, cosmological 
redshift, luminosity distance, single spin parameter, etc., and 
the uncertainties of corresponding parameters. The detection 
rates of our simulation correspond to a certain GRB population 
model. The population model and its hyper-parameters are 
followed from Yi et al. (2022) and Hendriks et al. (2023), 
which can successfully reproduce the detection rate consistent 
with historical GW and GRB observations (GWTC-3 and 
Fermi/Swift GRB catalogs). 

There is no capability currently for simulation of the arrival 
time delay between GW and sGRB signal by GWToolbox. 
Also, in our simulation, the sGRBs signal detected in joint 
detection with the GWs signal is on-axis. Although, observa- 
tionally, the Af bs of GW170817/GRB170817A is 1.7 s, some 
studies (like (Abbott et al. 2017b)) suggest that its corresp- 
onding sGRB is off-axis, and Afto»s is supposed to be smaller 
for on-axis cases. Therefore, in this work, we assume that the 
statistical fluctuation of the arrival time delay Afəbs is 
approximated by a truncated normal distribution between 
0.5s and 1.5s with a mean 1s, and a standard deviation 
of 0.5 s. 

Realization—After establishing a set of campaign para- 
meters such as the signal-to-noise ratio (S/N) threshold of GW 
detection, observation time, and GW detector, a catalog of 
samples is simulated. A set of Atop, is then randomly assigned, 
which is conveniently regarded as a “Realization” (i.e., 
simulating a realistic outcome once). Table 1 displays an 
example of a realization based on BNS merger events joint- 
detected by design aLIGO and Fermi/GBM over a 10 yr 
observation period with an S/N threshold of 8.0. Here we only 
show the important parameters used for calculation in the table. 
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Figure 1. The z distributions of the simulated catalog that are averaged over all realization, are presented in the top three panels for design aLIGO over a 10 yr 
observation period with S/N thresholds of 8.0, 7.0, and 6.0. The bottom three panels display the distribution for design aLIGO over a 20 yr observation period. 


The complete catalog also includes a single effective spin 
parameter, the uncertainty of masses of the binary stars, and so 
on. We can see that for each detected BNS merger event in 
Table 1, there is a redshift z and the corresponding uncertainty 
dz, as well as the arrival time delay At;5,. The uncertainties dz 
listed in Table 1 are estimated from those of the luminosity 
distance dL of the GW observation, which is calculated with a 
Fisher information matrix method with the GWToolbox (Yi 
et al. 2022). By applying the Monte Carlo (MC) algorithm to 
sample in the parameter space of z and Atops for each event, we 
obtain the distribution of the upper (lower) bounds. Both 
parameters are subject to normal distribution, with mean values 
taken from the catalog and standard deviations based on dz and 
0.01 s (which is considered as the uncertainty of the merger 
instant inferred from GW observation), respectively. Combin- 
ing all events in the catalog, we obtain the final bounds on óc/c 
using Equations (4) and (5) with this realization. Then we 
perform a new realization to obtain the new sample results of 
the catalog in the same observation setting to constrain the óc/ 
c. Instead of the Bayesian inference used by Nishizawa (2016), 
here we are using a simpler MC simulation, which allows us to 
effectively explore many realizations and different observa- 
tional campaigns. 


Taking into account the contingency of the distance of the 
detected events in a simulation, it is essential to conduct 
multiple realizations of a simulated observation setting. We call 
it a successful realization if there is at least one BNS merge 
event in which both the GW and the GRB are detected. After 
100 realizations are successfully performed, the results of all 
realizations are combined to obtain the final constraint interval 
of óc/c. 

The redshift distribution of samples that combines all 
realizations is shown in Figure 1, which uses design aLIGO 
over an observation period of 10/20 yr with different values of 
S/N threshold. Figure 2 illustrates the redshift distribution that 
combines all realizations from ET over three observation times 
of 3 months, 6 months, and 1 yr with S/N threshold — 8.0. As 
depicted in the distribution, nearly 2096 of the samples in 
design aLIGO observations have a redshift above 0.1, and most 
of the samples are greater than 0.5 for ET. Therefore, the 
measurement of the GW speed of all events takes into account 
cosmological effects. 


3. Result on óc/c 


The upper bound and the lower bound of 6c/c— the true 
intrinsic time difference range between GW and sGRB is 
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Figure 2. The z distributions of the simulated catalog that are averaged over all realization, which is observed by ET over different observation times with a S/N 


threshold of 8.0. 


Table 2 
The Constraints of Design aLIGO Joint-detected with Fermi/GBM in 10/ 
20 Yr Observation Time with Different S/N Thresholds 


Observation Time S/N Threshold Constraint of óc/c 


10 yr 8.0 (—1.88 x 1076, + 1.51 x 10717) 

10 yr 6.0 (—1.49 x 107!6, + 9.99 x 10718) 

20 yr 8.0 (—1.73 x 1076, + 1.31 x 10) 

20 yr 6.0 (—1.37 x 10716, + 9.40 x 10715) 
Table 3 


The Constraints of ET Joint-detected with Fermi/GBM in Different 
Observation Time with an S/N Threshold of 8.0 


Observation Time S/N Threshold Constraint of óc/c 


3 months 8.0 (—5.13 x 1077, + 1.27 x 10715) 
6 months 8.0 (—5.11 x 1077, + 1.15 x 10715) 
lyr 8.0 (—5.11 x 1077, + 1.08 x 10715) 


unknown. To obtain the upper bound of óc/c, we conserva- 
tively assume that the peak of the GW signal and the sGRB 
were emitted simultaneously, i.e., to consider the intrinsic time 
delay is 0, attributing the entire observed time delay to faster 
propagation of GW. For the lower bound, a conservative bound 
relative to the few seconds of intrinsic time delays has been 
discussed in (Nishizawa & Nakamura 2014; Li et al. 2016; 
Abbott et al. 2017b). We conservatively consider the intrinsic 
time delay to be 10s, i.e., the sGRB signal was emitted 10s 
after the GW signal. 

As discussed above, after performing one realization, we 
obtain a conservatively best constraint interval that combines 


^ The dispersion of the intergalactic medium (like plasma) has a negligible 
effect on the gamma-ray photon, with an expected propagation delay many 
orders of magnitude smaller than our assumption for Afin (Abbott et al. 2017b; 
Zhang 2023). 


all the events in the catalog. Combining all realizations then 
gives the final constraint range. To place a constraint on óc/c 
by combining all events joint-detected by design aLIGO and 
Fermi/GBM over an observation time of 20 yr, the best bound 
when S/N threshold is 6.0 is 9096 confidence interval 


-13705) x 10716 < a < +9.40+549 x 10-18. (6) 
For S/N threshold = 8.0, the number of events detected per 


realization decreases. We obtained the constraints (— 1.731337 x 


10-16, 1.317037 x 10717) of éc/c. Meanwhile, for ET, the bound 
over a time of 1 yr (S/N threshold = 8.0): 


IPM x 1071 < s £41080 x 1025. (7) 


The results above indicate that the constraint on óc/c has the 
potential to be about 10 times better than that of GW170817 
after 20 yr of observation by design aLIGO. Additionally, the 
constraint from ET in one year of observation time is only one 
order of magnitude stronger than the result calculated by design 
aLIGO in 20 yr, and two orders of magnitude stronger than 
GW170817. All of the constraints 90% confidence levels can 
be seen in Tables 2 and 3. 


4. Local Lorentz Violation 


SME is a comprehensive and effective eld theory that 
characterizes general Lorentz violations, which contains both 
general relativity and the Standard Model of particle physics. 
As a realistic theory, it can be applied to analyze observational 
and experimental data. The early foundational work can be 
found in (Colladay & Kostelecky 1998), while an updated 
review of observational and experimental results can be seen in 
(Kostelecky & Russell 2011). Additionally, some early study 
(Kostelecky & Mewes 2016; Kostelecky 2004; Bailey & 
Kostelecky 2006) showcases the foundational work in the 
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gravitational sector of the SME. The SME provides a general 
test framework for testing Lorentz invariance. One advantage 
of interpreting constraints on Lorentz invariance violation in 
terms of limits on the coefficients of the SME over model- 
independent tests is that it allows for direct comparison of 
results from different kinds of experiments (e.g., polarization 
and time-of-flight measurements) (Kislat & Krawczynski 2015). 
A Lorentz-violating term in the Lagrange density of the SME is 
an observer scalar density formed by contracting a Lorentz- 
violating operator with a coefficient that acts to govern the 
term. The operator can be characterized in part by its mass 
dimension d, which determines the dimensionality of the 
coefficient (Kostelecky & Mewes 2009). 

In a vacuum, Lorentz symmetry implies that all massless 
waves propagate at the speed of light. However, Lorentz 
symmetry is spontaneously violated when a medium is present, 
and propagation speeds can differ. Within a comprehensive 
effective eld theory description of Lorentz violation (Colladay 
& Kostelecky 1997, 1998; Kostelecky & Tasson 2015), there 
will be a Lorentz violation in both the gravitational and photon 
sectors of the SME. However, the two sectors of the violation 
are different, reflected in the effective-field theory operator by 
the fact that the coefficients of the Lorentz violation factor are 
different. The difference in the group velocity of GWs and EM 
waves is controlled by differences in coefficients for Lorentz 
violation in the gravitational sector and the photon sector at 
each mass dimension d. The coefficients for Lorentz violation 
can be described by an observer scalar and hence can be 
expanded in terms of spherical harmonics (Kostelecky & 
Mewes 2016; Kostelecky & Tasson 2015). At each fixed d, 
there are (d — 1) independent spherical coefficients, which are 
used to express the Lorentz-violating degrees of freedom. 

In the nondispersive, nonbirefringent limit at mass dimen- 
sion d — 4 (since every coefficient for Lorentz violation with 
d > 4 is associated with powers of momenta in the dispersion 
relation, only coefficients with d — 4 produce dispersion-free 
propagation (Kostelecky 2004)) using natural units and 
assuming that the nongravitational sectors, including the 
photon sector, are Lorentz invariant, attributing the difference 
in the group velocity to the Lorentz violation from the 
gravitational sector. In this case, the modified group velocity of 
v, takes the form (Kostelecký & Mewes 2016; Colladay & 
Kostelecky 1997, 1998; Kostelecky & Tasson 2015): 


y =14 ly cin, 3, (8) 


£m 


where Y,, is a basis of spherical harmonics, in which 
£<d—2=2, the spherical-basis coefficients 5 reflects the 
nine Lorentz-violating degrees of freedom for Lorentz violation 
in the gravitational sector. The direction (a, ô) refers to the sky 
position of the source. For the difference in group velocities, 
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we have 


Av = —» Yim(a, e (5c ors) 0) 


em, 


To place the constraints, we have proceeded under a 
simplified approach called a maximum-reach approach (Flow- 
ers et al. 2017) as used in Abbott et al. (2017b), separately 
considering each component in 5 and setting all other 
coefficients to be zero at one time. In the context of the 
maximum reach approach, sometimes referred to as a 
coefficient separation approach, a series of nine coefficients in 
Equation (9) are constrained one at a time using a single 
measurement in Abbott et al. (2017b). For the purpose of this 
demonstration of the principle, only the first isotropic item se 
is constrained in this paper, where the term is not associated 
with the sky direction of the event source. Calculating this 
leading term is sufficient for our purposes. The constraint of 3 
with different GW detectors (joint-detected with Fermi/GBM) 
are as follows (90% confidence level): 


design aLIGO: —1 x 10-55 < 5) <9 x 10-7, 
ET: —4 x 10-16 < $2 < 8 x 10-18. (10) 


These constraints are obtained from the observation period of 
20 yr (S/N threshold = 8.0) for design aLIGO and 1 yr (S/N 
threshold — 8.0) for ET. For design aLIGO, the isotropic upper 
bound of s(j shows a one-order of magnitude improvement 
compared to the constraint from GW170817 (Abbott et al. 
2017b), while for ET, it shows a three orders of magnitude 
improvement in only 1 yr observation period. 


5. Bound on Graviton Mass 


In theories of modified gravity (as described in (de Rham 
et al. 2017)), GWs generically propagate at a speed that differs 
from c and with frequency dependence dispersion relations. 
Will (1998) describes that in the case of inspiralling compact 
binaries, GWs emitted at low frequency early in the inspiral 
will travel slightly slower than those emitted at high frequency 
later, resulting in an offset in the relative arrival times at a 
detector. This modification affects the phase evolution of the 
observed inspiral gravitational waveform, similar to that caused 
by post-Newtonian corrections to quadrupole phasing. Matched 
filtering of the waveforms could limit such frequency- 
dependent variations in propagation speed and thereby 
constrain the graviton mass. 

If gravitation was propagated by a massive field, meaning 
the graviton is a massive particle, in which case the group 
velocity would be given by v=OE/Op, in a local inertial 
frame, 


2 
Vg g 

i= : 11 
E (11) 
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Figure 3. The frequency distribution of samples and the corresponding computed graviton mass distribution, from design aLIGO with S/N threshold equals 8.0 (top 


two panels) and 6.0 (bottom two panels) in the observation period of 20 yr. 


For the relativistic regime m « E, or in other words, if the 
frequency of the GWs satisfied Af > mg, where h is Planck 
constant, then the velocity of GWs (gravitons) will depend 
upon their frequency as shown in (Will 2014): 


2 
Ep tpm (12) 
c 2\Agf 


0.4 0.6 0.8 1.0 1.2 1.4 
m,/eV 1e-19 


where A, — h/m,c is the Compton wavelength of graviton. 
Therefore the bounds on speed of GW v, measured above are 
converted to the bound on graviton mass. 

In the frequency range centered around 100 Hz, the strain 
sensitivity of the design is ten times better than of the initial 
LIGO (Aasi et al. 2015). The signal of GW170817 detected by 
LIGO-Virgo Collaboration (LVC) has a frequency between 
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Figure 4. The frequency distribution of samples and the corresponding computed graviton mass m, distribution, from ET with S/N threshold equals 8.0 in the 


observation period of 1 yr. 


Table 4 
The Constraints on the óc/c, mg, and Soo from Different GW Detectors (Joint-detected with Fermi/GBM) with Varying Observation Periods and S/N Thresholds 


Detector Time S/N Threshold Constraint of óc/c Upper Bound of m,/eV S00 
Design LIGO 10 yr 8.0 (—1.88 x 10716, + 1.51 x 10777) 7.1 x 1077? (—1 x 107%, 1 x 1071$) 
10 yr 6.0 (—1.49 x 1076, + 9.99 x 10718) 6.2 x 10 7? (-1 x 1075, 7 x 10777) 
20 yr 8.0 (—1.73 x 10776, + 1.31 x 10777) 7. x 107? (-1 x 107 0,9 x 10777) 
20 yr 6.0 (—1.37 x 10776, + 9.40 x 10°18) 6.3 x 10 7? (—9 x 1076, 7 x 10777) 
ET 3 months 8.0 (—5.13 x 107", + 1.27 x 10715) 32 x 107? (—4 x 107!6, 9 x 10715) 
6 months 8.0 (—5.11 x 107", + 1.15 x 10718) 32 x 10 7? (—4 x 107!6, 8 x 10715) 
lyr 8.0 (—5.11 x 107", + 1.08 x 10°18) 32 x10 7? (—4 x 107!6, 8 x 10715) 


approximately 20 and 350 Hz (Abbott et al. 2017c), and the ET 
is designed to be 10 times more sensitive than advanced 
ground-based detectors, covering a frequency range of 1 to 10* 
Hz (Sathyaprakash et al. 2010). The limit on the GW speed 
corresponds to a narrow frequency range of GW in the merger 
stage. The corresponding frequency can be described by the 
following equation (Maggiore 2007): 


5/8 3/8 
fey (1) & 13an 12e) (3) (13) 


c T 


where 7 is the time to coalescence, and M, is the chirp mass of 
binaries. Equation (13) gives the GW frequency corresponding 
to r seconds before coalescence, with 7—0.01s as the 
uncertainty of the merger instant inferred from GW 


observation, or more specifically, the uncertainty of the peaking 
time of the GW signal. For each simulated GW source in the 
catalog, the corresponding binary mass data will be given. 
Substituting this data into Equation (13), we calculate the 
corresponding frequency for all samples from different 
observation campaigns. By then substituting the calculated 
frequency into Equation (12), we obtain the distribution of mz. 
Figures 3 and 4 show the frequency distribution and the 
corresponding computed graviton mass m, distribution of 
samples from design aLIGO and ET, respectively, which 
combine all realizations. 

At the 9096 confidence level, for design aLIGO in the 
observation period of 20yr, we obtain the constraint 
m, € 6.3 x 10 7^ eV. (S/N threshold = 6.0) and m, € 7.1 x 
10 ?? eV (S/N threshold = 8.0), respectively. For ET in the 
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observation period of 1yr, the bound of graviton mass is 
m,«32x10 7 eV (S/N threshold = 8.0), which yields 
almost the same constraints within the 10 ??^ orders of 
magnitude. 

In previous research, it is interesting to note that Nishizawa 
& Nakamura (2014) has placed a mass bound of the graviton 
within the 10 ?? order of magnitude in the absence of positive 
detection. This bound was derived from the study of 
observations of photons and neutrinos from compact binary 
mergers and supernovae. It is important to consider that 
Nishizawa & Nakamura (2014) uniformly chose the frequency 
of GWs to be around ~100 Hz for their calculations, while the 
frequency of GWs we calculated is almost located at ~700 Hz. 
This difference in frequency could indeed lead to variations in 
the constraints on the graviton mass. In addition, the graviton 
mass has been constrained by several observations of binary 
pulsars, resulting in a constraint of m, < 7.6 x 10 7? ev (Finn 
& Sutton 2002). Furthermore, Finn & Sutton (2002) has 
predicted that space-based GW detectors could place a better 
constraint on the graviton mass by observing multiple 
inspiralling black hole binaries, potentially leading to a 
constraint of up to m,~4x 10 ?6 ey (Berti et al. 2011). 
Overall, it is interesting to find that the bound on m, is 
comparable with that from binary pulsars, indicating the 
consistency and relevance of our results in the context of 
existing constraints. 


6. Conclusion and Discussion 


In this paper, a simulation study was conducted to constrain 
the GW propagation speed. The study focused on analyzing the 
difference in arrival time between GW and sGRB signals from 
compact binary mergers and considering the distance to the 
corresponding astrophysical event source. The simulation 
involved the joint observation of GW and sGRB signals from 
a BNS merger event using the GWToolbox. The study 
considered both second and third-generation ground-based GW 
detectors, specifically the aLIGO in design sensitivity and the 
ET, for joint detection with the Fermi/GBM in different 
observation settings. 

We have found that the joint detection of GWs and sGRB 
from the BNS merger has provided valuable insights into the 
constraining potential of aLIGO (in design sensitivity) and ET 
for the propagation speed of gravitational waves in future 
multi-messenger observations. Applying the MC algorithm, the 
upper bound of constraint óc/c reaches the order of 1077 for 
design aLIGO and 10 !? for ET, as indicated in Equations (6) 
and (7). Furthermore, the study obtained the coefficient 
—4 x 10776 < 5 — 8 x 107!8 for local Lorentz violation 
in the SME framework, along with a constraint on the graviton 
mass of mg <7.1(3.2) x 10 7? eV for design aLIGO (ET). 
Both of these values are constrained based on the GW speed 
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constraints from ET. Table 4 presents a summary of the results 
obtained under the various observation campaigns. 

The results show a significant improvement in the constraint 
of the GW speed, ranging from one order of magnitude for 
design aLIGO (over a 20 yr observation period) to three orders 
of magnitude for ET (over a 1 yr observation period), compared 
to the constraint of GW170817 in Abbott et al. (2017b). This 
highlights the advancements in the ability to constrain the GW 
speed through multi-messenger observations and the potential 
of advanced GW detectors, especially with the promise of 
future detectors such as ET. 

The bound on the GW propagation speed in this study is 
relatively conservative, as we simplified our analysis by using 
the same assumptions as in the calculation of Abbott et al. 
(2017b) for the upper and lower bound of óc/c. This involved 
assuming that the difference in intrinsic emission time of the 
two signals is 0s for the upper bound and 10s for the lower 
bound. By adopting this approach, we obtained the accuracy of 
this constraint that can be achieved in future observations. It's 
important to note that our simulations only considered the joint 
detection of GWs with on-axis sGRB from BNS merger events. 
However, as more joint detection events of GWs with other 
signals occur in the future, tighter bounds on the GW speed 
will be established, underscoring the importance of Multi- 
Messenger Astronomy in advancing our understanding of 
fundamental physics. 
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